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Abstract 

We develop a stochastic model to study the specific response of the immune system. The model is based on the 
dynamical interaction between Regulatory and Effector CD4+ T cells in the presence of Antigen Presenting Cells 
inside a lymphatic node. At a mean field level the model predicts the existence of different regimes where active, 
tolerant, or cyclic immune responses are possible. To study the model beyond mean field and to understand the 
specific responses of the immune system we use the Linear Noise Approximation and show that fluctuations due 
to finite size effects may strongly alter the mean field scenario. Moreover, it was found the existence of a certain 
characteristic frequency for the fluctuations. All the analytical predictions were compared with simulations using 
the Gillespie's algorithm. 
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1. Introduction 

The purpose of the immune system is to detect 
and neutralize the molecules, or cells, dangerous 
for the body without damaging the healthy cells. A 
fundamental process of the immune system is the 
maintenance of self-tolerance, i.e., the prevention 
of harmful immune responses against body compo- 



nents [1]. The biological significance of this process 
becomes very patent upon its failure during patho- 
logical conditions known as autoimmune diseases. 

The risk of autoimmunity cannot be dissociated 
from the capacity of the immune system to cope 
with diverse and fast evolving pathogens. The lat- 
ter is achieved by setting up a vast and diverse 
repertoire of antigen receptors expressed by lym- 
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phocytes, which as a whole is capable of recogniz- 
ing any possible antigen. Most lymphocytes have 
a unique antigen receptor (immunoglobulin in B- 
cells and TCR in T cells) that is encoded by a gene 
that results from somatic mutation and random as- 
sortment of gene segments in lymphocyte precur- 
sors. The randomness in the generation of antigen 
receptors makes it unavoidable that lymphocytes 
with receptors recognizing body antigens are also 
made. These autoreactive lymphocytes can poten- 
tially cause autoimmune diseases if their activation 
and clonal expansion are not prevented. The ques- 
tion is how is this avoided in healthy individuals? 

Regulatory CD4+ T cells, which express fork- 
head box protein 3 (FoxP3) are enriched in the 
CD25 pool of healthy individuals have been gain- 
ing increasing relevance in immunology [2]. Many 
lines of evidence indicate that these cells play a 
key role in the development of natural tolerance 
and in the prevention of autoimmune pathologies, 
by controlling the activation and proliferation of 
other autoreactive lymphocytes. The functional 
significance of these cells has broadened, as they 
were shown to modulate the immune response 
against pathogens, preventing the associated im- 
munopathology , and the rejection of transplants. 

Some years ago, Leon et al [3] proposed a simple 
model, consistent with these results able to explain 
the mechanism of immunological self-tolerance. 
The main hypothesis behind this model is that 
regulatory cells inhibit the proliferation of effector 
T cells, but depend on the latter for their own 
proliferation. The interaction between these cells 
is mediated by the presence of Antigen Presenting 
Cells (APC). 

The model, was developed to study the response 
of macroscopic populations of lymphocytes, i.e., the 
10 7 lymphocytes that may be present in a single 
lymph node and was technically approached study- 
ing the stability or numerically solving the corre- 
sponding differential equations [3] . This approach 
however, might be inappropriate to study the spe- 
cific clonal responses in the immune system. 

The 10 7 lymphocytes that exist in a lymph node, 
are really divided in approximately 10 5 or 10 6 dif- 
ferent clones [4]. I.e. sets of lymphocytes which rec- 
ognize the same antigen and participate in their 
own specific, rather independent, immunological 



response. Therefore, the number of cells involved in 
any particular clonal response is far smaller: 10 2 or 
10 4 T cells, which interact specifically with a small 
subset of all the antigen presenting cells available 
in the lymph node. One must naturally expect to 
found clones of different sizes, perhaps reflecting 
the ubiquitous nature of the antigen recognized by 
the T cells (i.e. the amount of APCs which are rec- 
ognized). However, the small size expected for most 
T cells clones would made differential equations an 
inappropriate tool to describe their dynamics. 

In this work, we first reformulated the model of 
Kalct and co. [3] within a more general stochastic 
frameworks of interactions between regulatory (R) 
and effector (E) T-cells. From this formulation we 
derive, on one hand a mean field description similar 
-but not identical- to the ones already known in 
the literature for the response of the entire immune 
system [5]. On the other, a self-consistent approach 
to unveil the role of the fluctuations in the specific 
response of the system. 

The rest of the paper is organized as follows. In 
the next section, that may be bypassed in a first 
reading of the manuscript, we present the main 
mathematical techniques involved in our work. 
Section 3 introduces and motivates the model. 
In section 4 we first compare, at the mean-field 
level, the predictions of our model with those of 
similar approaches in the literature. Then, we dis- 
cuss how fluctuations may affect these predictions 
highlighting their relevance for the clonal respon- 
seof the immune system. Finally, the conclusions 
and some remarks suggesting possible extensions 
of our work appear in Section 5. 



2. Mathematical Techniques 

When the effects of fluctuations are relevant, as 
it is often the case in biological systems, a modeling 
approach based on ordinary differential equations 
is not appropriate. 

This is well understood when the source of noise 
are the thermal fluctuations. But, while less men- 
tioned, it is also important when the noise is in- 
duced by the finite number of elements in the pop- 
ulations, i.e., molecules, cells, or individuals [6]. 
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A more convenient framework to describe these 
problems is a stochastic approach taking into ac- 
count the probabilistic nature of every process. The 
natural procedure is to write down a differential- 
difference equation that keeps track in time of the 
probability that the system would have a certain 
configuration (ft). Under the Markov assumption 
it is called the Master equation [7] 

J2 p (n,t)W,^, (f) 

where ft is a vector whose components are the num- 
bers of elements of every different species in the sys- 
tem, P(ft, t) is the probability of having the state 
ft at time t, and W~,^ n is the probability of the 
transition n' — > ft per unit time. 

Unfortunately, the Master equation can only 
be solved in special cases, usually, when the tran- 
sition probabilities are linear or constant in the 
state variables [8]. In more general situations we 
must, instead, resort to approximation methods. 
These approximation methods may be divided 
into two broad classes: numerical simulation algo- 
rithms, and perturbative calculations. Under the 
first class the Gillespie [9] algorithm has become 
the standard tool in the community (see appendix 
for a short description of the algorithm). It is a 
Monte Carlo's algorithm that simulates a single 
trajectory ft(t) consistent with the unknown prob- 
ability distribution P(ft,t) characterized by the 
Master equation. It guarantees to find the correct 
solution of the Master equation, but it is often 
time consuming and it is not suitable for param- 
eter exploration. On the contrary, perturbative 
methods lack the confidence given by a full and 
formal solution of the Master equation, but may 
provide important insight on the behavior of the 
system and on the relevance of the parameters. 

In this work we present two different perturba- 
tive solutions. The Linear Noise Approximation ( 
LNA) [10] which is an expansion around the large 
size solution of the corresponding Master equation. 
This gives us a tool to understand the relative size 
of the fluctuations and their relations with the pa- 



rameters of the model. On the other hand, the Ef- 
fective Stability Approximation ( ESA) [11], which 
is an extension of the LNA, provides a quantitative 
characterization of the effect of the fluctuations in 
the stability of the phases predicted by the corre- 
sponding mean field model. 

2.1. The Linear Noise Approximation 

An important difficulty in solving the Master 
equation arises from the discreteness of the vari- 
ables involved (ft). As long as the number of ele- 
ments increases the system evolution turns more 
regular and the mean-field equations give a more 
accurate description. The Linear Noise Approxi- 
mation is a systematic approximation method that 
rests upon the assumption that the deterministic 
evolution of the concentrations in the system can 
be meaningfully separated from the fluctuations, 
and that these fluctuations scale roughly as the 
square root of the system size fi. 

Moreover, in systems where the size of the pop- 
ulations differs in orders of magnitudes, it is ex- 
pected that the respective concentrations fluctuate 
within different scales. For this reason, we consider 
as many parameters Qj as species are involved in 
the reactions, and write the numbers of individu- 
als per species proportional to these values. Each 
fli will be a characteristic size scale for every par- 
ticular specie. 

Therefore, under the LNA, the population num- 
bers per specie are written as: 

ni = SliXi + (2) 

where x% is the deterministic prediction for the con- 
centration of the i th species with respect to the 
parameter f2 i; and an measures the fluctuations 
around x,. Note that formalizing the LNA in this 
way, there isn't any a priory assumption about the 
system size and concentrations higher than one are 
allowed. If for every specie 0, = Clj the standard 
formulation is recovered. 

Now we shall use the continuous variables on 
instead of the integers ni to write the probabil- 
ity distribution P(ft,t). Let us group the magni- 
tudes Xi into the macroscopic concentrations vec- 
tor x = (xi, x 2 , . . . , x;v), and the on into the fluc- 
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tuations vector a = (a\, . . . , ajv)- The matrix 
51 = diag(£li, Q2, ■ ■ ■ , &n) is a diagonal matrix 
whose principal elements are the extensive param- 
eters Cli with the units of the volume to which con- 
centrations X{ are refercd. 

Let us define the following variables: 



1 



uj = 



and ^ = JoioLi 



thus, we have 

rii = CtiXi + QiUjfli 



(3) 



(4) 



Taking (4) into account to expand the Master equa- 
tion (1) in powers of f^, at the leading order the 
following system of differential equations appears: 



dxj 
57 



M 

£ 



SijWj{Q.-x), 



(5) 



This corresponds to the deterministic rate equa- 
tions that are often used to describe the dynamics 
of such systems at a mean-field level. 

Here the time scale has been modified following 
t = Oit and Sij is the stoichiometry coefficient, 
the number in what the i th species varies when 
the j th process occurs. These coefficients will cor- 
respond with the element of row i and column j of 
the stoichiometry matrix S. 

In the next order approximation one finds: 



an(/3,^ir) 

a7 



N 



d\n0,si lT )p k ] 



+ 



i,k=l 



N 



where 



M 



A-ik — &i ^ ' Si 



i.k—1 



dWjin-x) 

dx k 



dpi 

d 2 U0, Pit) 

; dp t dp k 



and 



Dik = OiOk / SijSkjWj(Vt- : 



3=1 



(6) 



(7) 



(8) 



Aik and Dik are respectively the elements (i, k) of 
the matrices A and D, where, A can be identified 



with the Jacobian matrix of the system of differ- 
ential equations (5). 

Note that equation (6) is a Fokker-Plank equa- 
tion that characterizes the probability distribu- 
tion for the fluctuations ftir), centered on the 
macroscopic trajectory fi- x(t). 

The matrices A and D are independent of f3, 
which appears only linearly in the drift term. As a 
consequence, the distribution IT(/3, t) will be Gaus- 
sian for all time [12]. In particular, at equilibrium 
(x = x s ) the fluctuations are distributed with den- 
sity 



II s (/3) = [(2tt) n det[E]] 5 exp 



1 



7r 



E/3 



and variance S = ( (5 ■ (3 T ) determined by 



A-S + S-A T + D = 0, 



(9) 



(10) 



where A and D are evaluated on the studied fixed 
point [12]. 

The steady-state time correlation function is 



/3(r) ■ /F(t - t') ) = exp [At'} ■ E. 



(11) 



Moreover, if fluctuations are important it may 
be useful to study their properties. With this 
aim, a very powerful tool is the Fourier analysis 
of the equations that govern it. In our case it is 
the Fokker-Plank equation (6), who gives us all 
the information about the temporal behavior of 
fluctuations. Given the inconvenient form of this 
equation for our purpose, it is more reasonable 
to write it in a completely equivalent formulation 
more benevolent to investigation using Fourier 
transforms [13]. The problem can then be formu- 
lated as the set of stochastic differential equations 
of the Langevin type [12]: 



JV 



AikPk+Viir), 



(12) 



where /^(t) is a Gaussian noise with zero mean and 
correlation function given by 



(m(T)Vk(T')) = D lk S(r-T r ). 



(13) 



The power spectrum of the fluctuations can be 
found for every specie of the system by averaging 
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the square modulus of the Fourier transformation 
of In this way, 

Pi(u) = (|ft(a,)| 2 

N N 

= EE*« 1 M £> i*(* t )«M- ( 14 ) 

J=l fe=l 

Being $ lfe = zu;<S jfe - A ik and $^(cj) = $ ife (-w) 
[13] . With this expression we have an analytic mea- 
sure of the contribution of every frequency to the 
oscillatory behavior of the concentrations around 
each fixed point. Through it, the properties of the 
fluctuations can be studied. 

2.2. Effective Stability Approximation 

It is required just a bit of intuition to realize that 
the noise can affect the predictions of the deter- 
ministic rate equations models [11]. This is partic- 
ularly clear, once we question about the stability 
of the fixed point predicted by a mean-field model. 

In the deterministic analysis, the stability of the 
model ^ = f{x) to small perturbations is found 
by linearizing about the equilibrium point: x = 
x s + x p and the eigenvalues of the Jacobian J^ ' = 

|| _ ^ provide the decay rate of the exponential 

X — X B 

eigenmodes [14]. 

If we now consider the fluctuations around the 
small perturbation x p , the concentrations remain 



ill 

Ox 



x = x s + x p + w/3(t). The Jacobian matrix J = 
^ , will be expressed in terms of the 

=x s +u:/3(t) 

fluctuations lu(3(t) around the steady-state. In the 
limit lu — > 0, we can linearize J with respect to cu 
and the new Jacobian becomes 



T t I dJ 



= J (0) +wJ (1) (t), (15) 



Therefore, new effective eigenvalues A'i, includ- 
ing the influence of the intrinsic noise, characterize 
the stability of the phases [11]. 



A i — Aj + Ai 



(16) 



where Xi corr is proportional to oj, i.e., inversely pro- 
portional to the system size. 



When the system size is small, fluctuations are 
relevant, and in many cases the corrections Xi aorr 
can not be ignored. This small correction may, in- 
deed, change the sign of the eigenvalues (now A'j) 
and as a consequence, to change the stability of the 
phase. 



3. The cross-regulatory model 

In this section we motivate and present a 
stochastic model to describe the dynamics of two 
(possible small) populations of T cells E and R of 
the immune system. 

Some clues on how the regulatory CD4+ T cells 
suppress the response of other cells have been 
derived from well-correlated in vitro and in vivo 
experiments. These studies suggest that T-cell- 
mediated suppression is not mediated by soluble 
factors and require cell-to-cell contact [15] [16]. 
Moreover, regulatory CD4+ T cells can only sup- 
press the response of other cells if the ligands of 
both cells are expressed by the same Antigen Pre- 
senting Cell (APC). It also has been seen that 
regulatory T cells do not proliferate in vitro when 
cultivated alone in the presence of APCs. 

Despite these strong requirements, it is not yet 
clear what the nature of this mechanism is. The 
model requires a set of postulates that summarize 
the life cycle of these cells, and their interactions. 
Inspired by the success of Leon et al. [3] we are go- 
ing to assume that the proliferation of T cells oc- 
curs via its conjugation with an APC. An E (ef- 
fector) cell could duplicate only if conjugated with 
an APC which has no R (regulatory) cell. On the 
other hand, an R cell can duplicate only if con- 
jugated with an APC where at least one E cell is 
also conjugated. In that way, regulatory T cells ex- 
ist only in the presence of effector T cells, but the 
growth of the E cells is regulated by the presence 
of the R cells. A picture illustrative of these pos- 
tulates is shown in Figure 1. 

In vivo the total number of APCs and their ca- 
pacities to stimulate T cells may change, and it is 
perhaps a function of the T cells themselves, but we 
consider them as a constant population in which 
every APC has a finite and fixed number of conju- 
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Fig. 1. Cartoon illustrating the hypothetical mechanism 
presence-proliferation during simultaneous conjugation be- 
tween T lymphocytes and APCs. E cells promote R cells 
proliferation, while R cells avoid E proliferation. 

gation sites. Then, each site can be empty or occu- 
pied by a regulatory or effector T cell and all the 
cells that belong to the same phenotype (E or R) 
will have the same probability of being conjugated 
with any free site of any APC. 

Within these postulates we have to deal with 
the populations of free and conjugated cells for 
each specie, and the transition between the differ- 
ent states. This makes the model very cumbersome 
at the stochastic level, and almost intractable. For- 
tunately the processes of formation and dissocia- 
tion of the conjugates are relatively fast compared 
to the overall dynamics of the T cell populations. 
While a T cell remains conjugated with its APC 
for a few hours, the T cell mitotic cycle lasts on 
average 12 hours and the T cell lifespan is at least 
several days. Then, it is reasonable to assume that 
only free lymphocytes die, and that the numbers 
of conjugated T cells of both phenotypes are in 
quasi-steady state. 

Therefore, if a is the total number of APCs, and 
s, the number of conjugation sites per APC we may 
define: 



\ ~ E aB~ ) = Probabilit y P er unit time 
of combination of an E cell and an APC. 



k n E E ' 



k a f f E — = Probability per unit time of dissoci- 
ation of E to its conjugate. 

Here E and E b are the numbers of free and con- 
jugated E cells respectively, and k oriE and k ff B 



are positive parameters that characterize the pro- 
cesses of occurrence and disappearance of the con- 
jugates of this cellular species. 

In the same way, for R cells these probabilities 
per unit time are given by k onR R (^ as ~^ ls ~ R ~\ and 

k °ffn respectively. 



If now we set q e 



and q r 



the 



k °ff E , k °ffR ' 

stationary state for the formation and dissociation 
of conjugates for every species reads: 



E» 



asq e E 



1 + q e E + q r R 



R b 



asq r R 



l + q e E + q r R 



• (17) 



After this approximation, our analysis involves 
only four significant processes, the birth and death 
of effector and regulatory cells, and is far more 
tractable. The four processes are represented in 
Figure 2. 



Wi 



w 2 



w 3 



w 4 



APC 
. E cell 



Fig. 2. Reaction diagram indicating the events underlying 
the dynamics of APCs, E cells and R cells, as assumed in 
the model. 

Each event is characterized by an occurrence 
probability per unit time W(n — > fh) from the 
state n = (E, R) to the posterior one m = (E' , R'). 
It is, 



W(ft — > rh) 



= < 



Wi if fh= (E+1,R) 

W 2 if rn= (E,R + 1) 

W 3 if to = (E-1,R) 

Wi if to = (E,R-1) 



G 



and the stoichiometry matrix corresponding to this 
system, whose element contains the changes 
produced on the specie i th when it occurs the j th 
process is 



10-1 
1 0-1 



(18) 



Now, following the postulates of our model, a 
new E cell appears in the system if there already 
exists a combined E cell with an APC that has no 
R cells simultaneously conjugated. Then, the prob- 
ability of the increase of E in 1 is proportional to 
the number of APCs (a), times the probability to 
find in it a conjugated E cell (— ), times the proba- 
bility that this APC has no R cells simultaneously 
combined. The later has the form of an Hypergeo- 
mctric distribution [3] , but to gain in simplicity of 
the expressions we can manage it approximately as 

(l - , valid when a > 10 [17]. This takes 

us to the expression 



Wi = ip e a— 1 

as V as 



(19) 



The analysis for the growth of the R population 
in 1 unit is similar. A new regulatory cell appears in 
the system with certain probability if there exists 
already both R and E cells conjugated in the same 
APC during the same interval of time. The prob- 
ability of the existence of an E cell conjugated in 
a given APC where in addition exists space for an 

E" 



R cell can be written as 



1 



i 1 as ) 



Mul- 



tiplying this by the probability of randomly find 
that specific APC with an R cell combined in it 
gives the probability for this birth: 



W 2 = -0 r a— 
as 



as 



(20) 



In equations (19) and (20), ip e and ip r are param- 
eters for the effector and regulatory species that 
characterize the processes of formation and disso- 
ciation of the different conjugates. 

Finally, we assume that the death probabilities 
for each specie are proportional to the numbers 



of individuals, with proportionality factors m e for 
effector and m r for regulatory. 

Let us set e the concentration of E cells with re- 
spect to the quantity Q e , and r will be the concen- 
tration of R cells relative to fl r , where fi e and ft r 
are the characteristic mean values. The probabil- 
ity vector per unit time W, that contains the in- 
formation of the four processes of birth and death 
is formed by the elements: 



Wi = 

W 2 = 



1 + q e efl e 



^ e aq e eil e 

1 + q e efl e \1 + q e efl e 

ib r aq r rQ r 
. >< 

1 + q e efl e + q r rfl r 

1 + q r r£l r 



q r rQ, 



1 - 



1 + q e efl e + q r rQ, 7 



W 3 = m e e£L e 
Wi — m r rQ r . 



(21) 



(22) 

(23) 
(24) 



where the dependencies with E b and R b that ap- 
pear in (19) and (20) have been substituted by the 
expressions in (17). 



4. Results and Discussion 

4.1. Mean-field approximation 

Once we know all the elements of the transition 
probabilities per unit time, and the stoichiometry 
of the events, we are in conditions to write the 
deterministic rate equations that correspond to our 
model (see eq. (5) ). 



dr 



dr n e n e 
d-T = n w >-n w ^ 



(25) 
(26) 



where the time has been measured according to 
t = Q e r. 

In the steady-state, ^ = ^ = 0, the solutions 
for the deterministic system are given by a fixed 
pair (e*, r*) of numbers. Depending on the param- 
eters, the stationary solutions with biological sense 
may be of three types: (0, 0), (e* > 0, 0) and (e* > 
0,r* > 0). 



7 



Attending to the stability of these solutions we 
define five phases that are sketched in Figure 3. 



Stability according to eigenvalues 



12 
10 

8 
6 
4 








0.7 1.4 2.1 2.8 

Vr 



3.5 



Fig. 3. Keeping constant the parameters s = 3, a = 1000, 
n e = 1000, Q r = 2000, q e = 10~ 3 , q r = 5 ■ 10~ 4 , m e = 1.0 
and m r = 0.1, the five parameter regions with different 
kinds of stability depending on the [^j r ;i/j e ] pairs arc shown. 
Each ordered pair was classified attending to the signs 
of the eigenvalues of the Jacobian evaluated on the fixed 
points. Solid lines represent the region borders. 

Phase 0: In this phase only the trivial fixed 
point (0, 0) exists. It is characterized by the rela- 
tion tp e ciq e < m e and within this phase, this fixed 
point is always stable. 

Phase I: In this phase, only the point (e* > 0, 0) 
is stable. 

Phase II: Only the point (e* > 0, r* > 0) is 
stable. 

Phase III: It is a bistable phase. Both points 
(e* > 0, 0) and (e* > 0, r* > 0) are stable 

Phase IV: In this phase, the fixed points pre- 
sented above are unstable. A limit cycle develops. 

This phase diagram is very similar to the one 
reported in [3]. However, for completeness and be- 
cause the Phase IV was not predicted before we 
discuss in detail below the biological interpretation 
of the five regions. 

Since E cells themselves do not mediate the ef- 
fective immune response, but they trigger it, we 
interpret the state dominated by regulatory T cells 
(e*,r*) as tolerance or unresponsiveness, and the 
state dominated by effector (e* > 0, 0) as an effec- 
tive immune response or autoimmunity. Such in- 



terpretation arises from the fact that when R cells 
exist in the system, this population competes with 
the effectors for the APCs, limiting the growth of 
the E population. 

Particularly relevant from the immunological 
point of view is Phase III. Operating in this re- 
gion, the model can switch between the states 
dominated by regulatory or effector T cells in a 
reversible way [3] . 

For the fixed points of the kind (e* > 0, 0), e* = 
^m^Ta"" ' the eigenvalues of the Jacobian matrix 
are: 



Al(e*,0) = m e^e + 



aq e ip e 



(27) 



and 



^2( e *,0) — — O e m r + 



Q e q r i[) r 



q e ip e \aq e ip e 



(28) 



the first one of these eigenvalues is always nega- 
tive and its corresponding eigenvector is Vi( e *,o) = 
( 1 , 0) . It guarantees that , if not perturbed along the 
r direction, the fixed point is stable. This means 
that in those phases where both species can coexist 
(77, 777 and IV), once the R population is extin- 
guished, to take the system away from the (e*,0) 
state it is necessary that at least one R cell comes 
from outside the system. The absence of such a 
perturbation can be interpreted as the lost of tol- 
erance, and the appearance of an effective immune 
response. It explains the fact that thymectomized 
animals are more susceptible to procedures that 
induce autoimmunity that euthymic animals [18]. 

It is also well known that if in a tolerant organ- 
ism the number of APCs is increased enough, the 
immune system may turn on an immunological re- 
sponse. To study under which conditions our model 
reproduces this effect it is useful to understand 
that to change the parameters space from a state 
given by (a,tp r ,tp e ) to another characterized by 
(na, ip r , ipe) is equivalent, in terms of the determin- 
istic equations, to change the state (a,ip r ,tp e ) to 
(a,nip r ,nip e ) along the straight line that matches 
them. So, the increase in the number of APCs (a) 
may be interpreted as a motion along a straight 



8 



line with slope ^jf- that passes through the origin 
in the ip r vs ip e map. 

In our model, a tolerant organism is located in 
phases II or 777. In the case of a system that starts 
in region III, raising the number of APCs makes 
the system first, to evolve into region I (see the 
line with points in Figure 3). Biologically, it may 
be interpreted as the entrance into a region where 
the tolerance disappears and the animals develops 
immunological response. On the other hand, when 
the organism starts in the phase II, increasing the 
number of APCs drifts the system into a cyclic 
immune response (see the dashed line in Figure 3). 

It is worth to highlight that our mean-field pre- 
dictions compare very well with those obtained by 
K. Leon et al.[3] for a similar model of R-E in- 
teractions regulated by APCs. However, our phase 
diagram is a bit richer. Exploring a parameter re- 
gion where q e ^ q r the bi-stable phase III may be 
now bounded by a phase IV where the biologically 
meaningful fixed points are unstable. Using param- 
eters in this zone we have numerically solved the 
system of differential equations (25) and (26). As 
Figure 4 shows, in this zone appear very well de- 
fined closed orbits, representative of limit cycles. 

Numerical solution in the limit cycle 




R cells concentration 

Fig. 4. With the parameters of Figure 3, and the pair 
[ipr = 2.2; ip e = 11-0], it is shown the numerical solution for 
the deterministic rate equations (25) and (26). The trajec- 
tory surrounds the unstable fixed point (e* , r*)(represented 
by a white circle) and evolves to a cyclic behavior. 

The existence of this limit cycle allows the model 
to reproduce the emergence of a cyclic immune re- 



sponse in pathologies such as Multiple Sclerosis 
[19] , [20] . In this way, changing the parameters one 
can observe two different kind of immune response, 
one stationary and another cyclic. Both patologies 
are observed in the Multiple Sclerosis, perhaps be- 
cause genetic variations between individuals turn 
on the immune response within a different set of 
parameters. To our knowledge, no other mathe- 
matical model based on the dynamics of T cells 
regulated by R cells had predicted it together with 
effective response and tolerance. 

Moreover, it is useful to compare the previous 
mean-field results with our original stochastic 
model. The reasons are two-fold, first it should 
confirm that the mean-field approach is correct 
when studying large systems and second, it should 
shed some light on the behavior of smaller systems, 
like those characteristics of specific responses. To 
do this, we run Gillespie's algorithm. In Figure 5 
we show the results of typical runs for parameter 
values near the different fixed points, while in Fig- 
ure 6 the numerical solutions of the deterministic 
equations for e(r) and r(r) are compared with 
the Gillespie's simulations for the cyclic regime. 
These graphs demonstrate that, within a stochas- 
tic approach, strong fluctuations appear around 
the values predicted by the mean-field equations. 
Below, we show that if the size of the population 
is small enough these fluctuations become relevant 
and may qualitatively change the results of the 
mean field approximation. 

4.2. The role of the fluctuations 

In smaller systems the noise may induce qualita- 
tive changes in the mean field dynamics discussed 
above. This fact might be particularly relevant, 
since it could imply that the dynamics of small T 
cells clones would be significantly different from 
the one described here and in previous works [3]. 
We have carried out a series of Gillespies simula- 
tions with our stochastic model using parameters 
sets that belong to region III when classified ac- 
cording to the mean-field results. Our results illus- 
trated that steady state stabilities, specially for the 
steady state dominated by regulatory T cells (tol- 
erant steady-state), could be significantly affected 
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Parameter region Parameter region I Parameter region II 




Fig. 5. Keeping constant s = 3, a = 1000, Q e = 1000, Q r = 2000, q e = 10~ J , q r = 5 • 10 , m e = m r = 1, the figure shows 
results of simulations using the Gillespie's algorithm, for pairs of parameters [ip r ;tf> e ] corresponding to different regions of 
the phase diaagram. Panel a:[i/i r = 185.6 ; ip e = 78.0], in region O, the extinction of both populations is verified. Panel 
b:[i/v = 4.0 ; ip e = 110.0] in region /, just the effector population survives. Panel c:[ip r = 185.6 ; ip e = 198.0] in region II, 
both species, effector and regulatory, persist. Solid lines represent the deterministic values for fixed points, around which 
concentrations fluctuate. 



Numerical solution in the limit cycle 
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Fig. 6. With the same parameters of Figure 3 and the pair [ip r = 2.2; %f> e = 11.0] from region IV, the numerical solution 
for the deterministic rate equations (25) and (26) is compared with a simulation made with the Gillespie's algorithm. The 
trajectory surrounds the unstable fixed point (e*,r*) and evolves to a cyclic behavior. Both the amplitude and period of 
cycles fluctuate in the simulations. 



by stochastic fluctuations. Figure 7 shows how in 
simulations the T cell concentrations remain oscil- 
lating around the mean-field value of the tolerant 
steady-state where the simulations began. However 
after a while fluctuations drive the system away 
into a state dominated by effector T cells (immune 
steady state). In other cases simulations shows that 
the tolerant steady state became inestable under 
stochastic fluctuations but it dynamically evolve 
into a trivial steady state with zero effector and 
regulatory T cells. 



A new phase diagram (predicted with the help 
of (16) is represented in Figure 8, and has some dif- 
ferences when compared with Figure 3. The more 
striking feature of this map is the apparition of a 
new zone embedded in region IV with the charac- 
teristics of phase II. In this region of the parame- 
ters, instead of having a cyclic behavior, the system 
evolves into a tolerant state (region II) . This con- 
clusion can be corroborated by simulations, and is 
in contradiction with what the normal analysis of 
eigenvalues predicts. 
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Fig. 7. With the fixed parameters of Figures 3 and 8, is 
shown a simulation using the pair of parameters [ip r = 1.17 
; ip e = 8.29] that belong to region 777 according to the nor- 
mally calculated eigenvalues (Figure 3), and belong to re- 
gion I according to the corrected ones (Figure 8). The sim- 
ulation starts on the corresponding fixed point (e* , r * ) that 
was supposed to be stable, but fluctuations have changed 
its stability. 

This result sustains the idea that, increasing the 
number of APCs, a cyclic immune response may 
be turned into a tolerant response (see the dashed 
line joining regions IV and II in Figure 8). Along 
this line, when the number of APCs is increased, 
simulations show a stretching of the amplitude of 
the cycles until they reduce to a mere fluctuating 
behavior around the deterministic fixed point as in 
Figure 5 c. 

On the other hand, fluctuations, also reduced the 
bistable region III favoring the immune response 
of the system. Note however, how small spots of 
bistability persist within region I. 

Moreover, we study the effect of these fluctua- 
tions when the system is subject to external per- 
turbations. For example, in the phase III, the sys- 
tem, depending on the initial conditions may stay 
in the tolerant state (e > 0, r > 0) or in the state 
(e > 0, r = 0). In the mean field model, to switch 
the system from one state to the other, it must be 
strongly perturbed and moved away from the re- 
spective basin of attraction of the fixed point. In 
the case of small systems, small perturbations may 
switch the state of the system, provided they are 
frequent enough. 
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Fig. 8. With the same parameters of Figure 3 to build this 
new map, there were calculated for each point [i/y;i/> e ] the 
corrected eigenvalues A'j = A^ + ^i corr , and according to 
the signs of these, it was classified every point in one of the 
regions O to IV. Black points correspond to region 777. 
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Fig. 9. Transition from the tolerant state (e > 0, r > 0) 
to the responsive state (e > 0, r = 0) as a function of the 
number of external perturbations when the system is in 
region III . 

In figure 9 we show the probability (P c ) of 
switching from the tolerant state to the state 
(e > 0, r = 0) as a function of the number of 
perturbations. In this case, with a period T the 
number of cells at a given time was reduced an 80 
percent. T s is the size of the simulations (different 
symbols mean three different T s ). As can be easily 
see if the number of perturbations is small enough 
T s /T — > 0, the system do not switch from one 
state to the other. It is always tolerant. On the 
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other hand, when T s /T — > oo, i.e. when the period 
of the perturbations is small enough, the system 
goes to the responsive state with probability 1. 
The existence of a clear threshold for the mini- 
mum number of perturbations needed to change 
the state also signatures the importance of the 
fluctuations for the ocurrence of this transition. 

All these results confirm our previous intu- 
ition. The role of the fluctuations may be rele- 
vant and may alter significantly the predictions 
of the mean-field models. In particular, regions 
that were considered bi-stable may now develop 
immune response, and regions in which a cyclic 
pattern was predicted may become fully tolerant. 
Moreover, also the response of the system to ex- 
ternal perturbations change. Perturbations, that 
may be produce only transient behaviours in the 
deterministic model, may become relevant when 
fluctuations are taken into account. Of course, the 
magnitude of these results are parameter depen- 
dent and probably richer diagrams may be found 
exploring further the model. But our results al- 
ready prove that the interpretation of the specific 
immune response in base of mean field models 
should be done with care. 



4.3. Power spectrum 

Any single run of the stochastic system (as pan- 
els b and c in Figure 5) shows that coherent oscilla- 
tions are sustained in parameter regimes in which 
the deterministic equations approach a fixed point. 
It evidences the fact that fluctuations can not be 
safely ignored in systems composed of a few thou- 
sands of elements. 

Therefore, it is important to understand, not 
only the consequences of these fluctuations as done 
above, but also their properties. To do this, we ana- 
lyze the characteristics of these fluctuations around 
the steady state, using equation (14): 

N N 
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Then, evaluating A and D in the corresponding 
fixed point, and using them to calculate eq. (29) we 
can determine the contribution of every frequency 
to the composition of the oscillations. Then, the 
discrete Fourier transform of a data of concentra- 
tion and time that comes from the simulations can 
be compared with the power spectra curve calcu- 
lated with (29). It would give us a way to know 
how good the made approximations were. 

In Figure 10 we show the power spectra calcu- 
lated by (29) and the one obtained by averaging 
the Fourier transforms of 6000 different Gillespie's 
simulations. Both have been normalized in a way 
that P(lj = 0) = 1 and it is evident that the coin- 
cidence is remarkably good. 

The peak in the power spectrum indicates the 
existence of a resonance frequency. In it, the fluc- 
tuations of the concentrations that define the state 
of the system amplify their values as a function of 
the parameters of the model. This effect may be 
interpreted as the induction by the internal noise 
of a a kind of stochastic resonance that is respon- 
sible for the emerging of amplified oscillations in 
systems which, when described through the deter- 
ministic rate equations for a given set of parame- 
ters lack a cyclic behavior. 

From the technical point of view, this character- 
istic frequency arises because of the existence of 
complex eigenvalues with non-zero imaginary parts 
for the Jacobian at the fixed point. In the stan- 
dard stability analysis, the imaginary part of the 
eigenvalues are responsible for oscillatory transient 
behaviors near the fixed points. However, here the 
stochasticity due to the smallness of the system 
results in persistent perturbations away from this 
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Fig. 10. Superposition of the power spectrum for E cells 
calculated with the analytic expression (29) and the square 
of the absolute value of the Fourier transform of 6000 
Gillespie's simulations. The normalization condition is 
P(ui = 0) = 1. The used parameters are [ip r = 185.6 
; Ve = 198.0]; q e = 10~ 5 ; q r = 5 ■ 10~ 4 ; Hi = 1000; 
Q2 = 2000; m e = m r = 1; s = 3 and a = 1000, corre- 
sponding to region 77. 

fixed point, so that both features together result 
in an overall oscillatory effect. 

The existence of this characteristic frequency 
could be useful to induce transition between two 
states of the system just by making a parameter 
(such like the number of APCs) to fluctuate with 
a small amplitude but with the appropriate fre- 
quency. Simulations corroborate this when it is in- 
duced a change from a tolerant phase to an immune 
one by making a to fluctuate near a III-I region 
border. Work is in progress to check the reliability 
of this approach to develop new therapeutic tech- 
niques. 



5. Conclusions 

We present a stochastic model of interactions 
between Effector and Regulatory CD4+ T cells in 
the presence of Antigen Presenting Cells that is 
able to show active, tolerance, or cyclic immune re- 
sponses. We characterize the phase diagram in the 
mean field approximation, and then analyze the 
corrections, due to finite size effects, to this phase 



diagram in the presence of fluctuations. This cor- 
rected phase diagram may be understood as a first 
attempt to characterize a stochastic model for the 
specific response of the immune system. In fact, we 
prove that the fluctuations may strongly alter the 
mean field predictions turning for example toler- 
ant zones in responsive or cyclic responses into tol- 
erant, and therefore that any analysis of the clonal 
response of the immune system based on mean field 
predictions must be taken with caution. Moreover, 
we show that our model present sustained oscil- 
lations at a characteristic frequency that may be 
relevant to understand, or treat these clonal re- 
sponses. 



Appendix A: Gillespie's Algorithm 

The occurrence of every process that can take 
place into our system has a very strong stochas- 
tic component. Every reaction depends on a lot of 
variables that can only be described in a probabilis- 
tic way, so in a given state, both the next process 
that would take place, and the moment in which 
it will occur are random variables. Even when it 
is not possible to find an analytic solution for the 
Master equation that governs the system, micro- 
scopic simulations of the processes defined by the 
relations (21) to (24) can be carried out using the 
algorithm originally proposed by Gillespie [9] . 

The algorithm keeps the random fashion in 
which processes occur. The backbone of such a 
simulation consists in determining which one will 
be the next process, and when will it occur. Every 
run of the simulation will give a different temporal 
sequence of concentration that is, according with 
the Master equation of the system. 

The procedure can be synthesized in three main 
steps: 

(i) Calculate the occurrence probabilities (Wj) 
per unit time for every process that can take 
place in the system. 

(ii) Generate two random numbers n and r 2 us- 
ing a unit-interval uniform random number 
generator, and calculate At and \i according 
to (30) and (31): 
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At = 4t In | — ) and (30) 
W \ri ) 

H-l fi 

Y,W,<r 2 W <J2 W ^ ( 31 ) 
where 

M 

W = J2Wj. (32) 

(iii) Using the At and £i values obtained, increase 
the time t by At, and adjust the cellular pop- 
ulation levels to reflect the occurrence of the 
[i th process. 
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